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Abstract. We consider the problem of optimizing the locations of source nodes in 
transportation networks. A reduction of the fraction of surplus nodes induces a glassy 
transition. In contrast to most constraint satisfaction problems involving discrete 
variables, our problem involves continuous variables which lead to cavity fields in the 
form of functions. The one-step replica symmetry breaking (IRSB) solution involves 
solving a stable distribution of functionals, which is in general infeasible. In this 
paper, we obtain small closed sets of functional cavity fields and demonstrate how 
functional recursions are converted to simple recursions of probabilities, which make 
the IRSB solution feasible. The physical results in the replica symmetric (RS) and the 
IRSB frameworks are thus derived and the stability of the RS and IRSB solutions are 
examined. 
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1. Introduction 

Constraint satisfaction problems (CSPs), which are highly relevant to many applications 
such as electronic circuit design and frequency assignment in cellular mobile networks, 
have been studied in the fields of applied mathematics, computer science and 
engineering. Despite their usefulness, many CSPs are NP-complete problems [1] 
associated with algorithmic hardness. It is thus important to understand the physical 
origin of their hardness and map out the easy and hard regimes for typical instances. 
Physicists approach the problem by making analogy between CSPs and spin glasses 
[21 [3]. Objective functions are mapped to spin glass Hamiltonians, enabling the analysis 
of CSPs using statistical physical techniques. Successful examples are found in the K- 
satisfiability problem [1] , graph coloring [5] , and vertex cover [Gj [TJ [8] . They suggest a 
rich physical picture of CSPs corresponding to the glassy phase in spin glasses. 

In this paper we study a problem having a wide range of applications and sharing 
the characteristics of CSPs. Specifically, we consider the optimal locations of source 
nodes in transportation networks. Transportation networks consist of nodes with either 
surplus or deficiency of resources, and an important problem is to distribute them so as 
to achieve a networkwide satisfaction with a minimum transportation cost [HI [ini IH] ■ 
This problem is important in load balancing in computer networks [12] and network 
flow of commodities [13]. Progress has been made in generalizing the message-passing 
technique of discrete variables to the passing of cavity energy functions in terms of the 
continuous current variables [HI [THl E] ■ 

Subsequent work considered networks in which shortages are allowed but cost 
penalty is imposed [I4j. This models applications such as communications networks 
where shortages are detrimental to the performance of the nodes. Their effects were 
modeled by step-like shortage costs. This high nonlinearity gives rise to unique behavior 
and a physical picture absent in the previous models. When the shortage cost is 
comparable to the transportation cost, the total cost may be optimized either by saving 
the transportation cost feeding a poor node while sacrificing the satisfaction of the 
node, or by saving the shortage cost while spending more on the transportation cost. 
The picture in reminiscent of the learning of noisy examples in perceptrons, where the 
field distribution of the examples consist of the bands, corresponding to the learned 
and sacrificed examples respectively [151 |T6l [T71 [H] . As a result, frustration arises from 
competition for resources among connected nodes. Numerous metastable states emerge, 
leading to typical glassy behavior. 

The problem of optimal source location in this paper addresses an even more general 
and practically relevant issue in network design and optimization. Compared with [H] 
where some nodes remain unsatisfied in the optimized state, this paper moves one step 
forward and considers the situation in which the location of the source nodes can also be 
optimized, and all nodes are satisfied. The source location problem has wide applications 
in the design of optimized transportation networks. For example, the optimal locations 
of access points in wireless networks can be determined by balancing the signaling cost 
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of the access points and the power and bandwidth hmitations of the channels hnking 
the mobile subscribers (which can be expressed as the transportation cost). 

As demonstrated in P, [101 [H], the resource allocation problem involves passing 
messages of continuous variables. When the cost function includes nonlinear terms, the 
messages generally become extremely complicated. However, as will be described in 
this paper, there are phases where the space of continuous messages can be replaced by 
small closed sets of cavity energy functions, and their recursions can be converted to 
simple recursions of probabilities. In the context of the source location problem, this 
takes place when the consumer nodes form small clusters surrounded by source nodes. 
When the ratio of the installation cost of the source nodes and the transportation cost 
changes, regimes with different maximum cluster sizes are observed, resembling the 
Devil's staircase observed in the circle map and other dynamical systems [26]. We will 
show that the use of small closed sets of cavity energy functions is particularly successful 
in the singlet regime where the consumer nodes are isolated (clusters of size 1) and, when 
the cost ratio is commensurate, in the doublet regime where the consumer nodes can be 
paired or isolated. 

In fact, clusters formed by similar energetic considerations have been found to play 
an important role in disordered systems such as the random field Ising model (RFIM) 
|19j . as illustrated in Fig. [1] Indeed, domain sizes in RFIM are determined by the 
interplay between the random field energy and the domain wall energy, giving rise to the 
so-called Griffiths singularities and cascades of phase transitions [201 [21]. Analogously, 
cluster sizes in the source location problem are determined by the balanced between the 
installation cost of the source nodes and the transportation costs. 

Another class of problems that exhibits similar cluster behaviors is the lattice glass 
models [221 [23], |2l] . In the lattice glass models, each site of a network can be occupied or 
empty, but the number of nearest neighbors cannot exceed a maximum. Analogously, 
the energetics of cluster formation in the source location problem implies that two 
neighboring consumer nodes have effective repulsions, since their contiguity prevents 
them from drawing resources from more neighboring source nodes. However, the present 
model is richer in behavior, since the cluster energy depends on the current distribution 
in non-trivial ways, rather than merely counting the number of neighbors. When the 
particle density in a lattice glass increases, configurations of clusters are formed, causing 
the system to undergo dynamical freezing transitions preceding the ideal glass transition. 
Similar transitions will be reported in the source location problem. 

The vertex cover problem [6l [8] , one of the NP-complete problems in computational 
complexity theory [1] that attracted recent attention from physicists, also shares similar 
cluster behaviors. Drawing the analogy of assigning guards in a museum [6], each site 
of a network can be covered or uncovered, but none of the links can have both vertices 
uncovered. Hence the uncovered nodes form a configuration with an effective repulsion 
among them. Indeed, an equivalent configuration of consumer nodes can be found in 
the singlet regime of the source location problem, since for neighboring consumer nodes, 
each has one less link to draw resources from the neighboring source nodes. When the 
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Figure 1. (a)-(b) RFIM of two spins and four fixed boundary spins. Symbols: filled 
(unfilled) circles -up (down) spins, arrows - direction of random fields of magnitude h, 
thick segments - frustrated couplings of magnitude J. When h/J increases from below 
2 to above, pattern (a) changes to (b) to form a cluster demarcated by a domain wall of 
frustrated couplings, (c)-(d) The source location problem defined by Eq. ^ with two 
nodes and four fixed boundary nodes, and — —1 for all i. Symbols: filled (unfilled) 
circles - consumer (source) nodes, arrows - current flows. When decreases from 
above to below, pattern (c) changes to (d) to form a cluster demarcated by a 
"domain wall" of inward current flow. 



average connectivity of the network increases, we will show that it undergoes a glassy 
transition analogous to that in the vertex cover problem. 

The origin of these interesting phenomena can be traced to the presence of 
frustrations, which refer to the conflicts between competing interaction energies in the 
system [25]. This connects our problem with a broad class of network CSPs in which 
frustrations are inherent. When the system consists of numerous states, the replica 
symmetry-breaking (RSB) solution is applicable, but the distribution of the cavity 
energy functions is in general infeasible to solve. Nevertheless, with the small closed set 
of cavity energy functions introduced in this paper, the one-step RSB (IRSB) solution 
becomes computationally feasible. The physical results based on the replica symmetric 
(RS) ansatz and the IRSB configurational entropy are thus derived. 

The paper is organized as follows. We introduce our problem in Section [21 
presenting simulation results of the optimal behaviors. In Section [3l we discuss the 
general RS formalism and derive the piecewise quadratic ansatz of the cavity energy 
functions. In Section [H we obtain a small closed set of cavity energy functions in the 
singlet regime and demonstrate how functional recursions are converted to recursions 
of probabilities. The average energy, the fraction of soft nodes and the RS-RSB phase 
diagram are derived. In Section [6], we apply the small closed set of cavity fields to 
the IRSB formalism and obtain results of the configurational entropy. The conclusion 
is given in Section [71 In Appendix A we describe how the small closed set of cavity 
energy functions can be applied to the doublet regime in the commensurate case, and 
the corresponding RS and IRSB results are derived. 
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A,, 


before optimization 




after optimization 


> 


surplus node 


> 


source node 
(resource 
providing) 


< 


deficient node 




> 


consumer node 
(consuming) 



^The deficient nodes with < after optimization are converted to source nodes. 
Table 1. A summary of glossary used in this paper. 

2. The Model 

2.1. Model Formulation 

We consider a network of nodes, labelled i = 1 . . . N. Each node i is connected 
randomly to a set Ci of K neighbors. Each node i has capacity Af, nodes with 
positive and negative values of Aj correspond to surplus and deficient nodes respectively. 
The capacities Aj are randomly drawn from a distribution of p(Aj). With network 
applications in mind, we consider a bimodal distribution in which Aj = y4(^ 1) with 
probability 0s > and Aj = — 1 with probability (j)^ = I — 4>s- Naturally, the surplus 
nodes serve as source nodes providing resources to the consumer nodes. However, to 
minimize cost functions that include transportation costs, it is often desirable to convert 
some deficient nodes into source nodes as well. Hence in general, the task is to optimally 
locate these extra source nodes so as to minimize the total cost function. The relevant 
glossary used in this paper is summarized in Table [H 

We first consider the minimization of the cost function in [T3], whose optimization 
variables are the currents i/ij = —i/ji of real values from node j to node i, 

C,i = Ai + J^jedyij final resource of node i, and Q{x) = 1 when a; > 0, and 

otherwise. The link connecting nodes i and j is denoted as (ij). The first term 
corresponds to the unsatisfaction cost imposed on nodes with negative final resource. 
The second term is the transportation cost. This cost function models load balancing 
situations in which insufficient provision of resources to a deficient node produces 
detrimental effects on it (irrespective of the magnitude of insufficiency). 

The key to applying the cost function in Eq. ([1]) to optimize the location of 
source nodes is to note that once the final resource of a deficient node is negative. 
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the unsatisfaction cost remains the same even when its resources are maximally drawn 
by other nodes of the network. Hence the deficient node effectively becomes a resource 
provider. Ref. [H] contains many such examples. If we consider the coefficient to 
be the installation cost of a source node, then we can solve the optimal source location 
problem by first minimizing the cost function in Eq. ([T]), then identifying the deficient 
nodes whose final resources are negative, and converting them to source nodes. 

Formally, in the optimal source location problem, we introduce the state variables 
Sj = ±1 for deficient nodes when node i is a consumer or a source node respectively. 
The cost function is then 



subject to > for Si = 1. No constraints are imposed on nodes with Si = —1, since an 
arbitrary amount of resource can be provided when they are converted to source nodes. 
Afo is the set of deficient nodes. 

To check the equivalence between the cost functions in Eqs. (II]) and ([2]), we can 
easily see that when Sj = 1, the installation (or unsatisfaction) cost vanishes in both cost 
functions. When Si = —1, we only have to consider the case < 0, and the installation 
(or unsatisfaction) cost is in both cost functions. This is because when > 0, we 
can set Sj = 1 to minimize the total cost. 

Note that the cost function of the optimal source location problem is identical to 
that in [14J, but the interpretation is far more relevant to network applications. All 
previous results on networks with nodes of negative capacity can be directly mapped to 
networks whose unsatisfied nodes are replaced by source nodes. For example, the single- 
sat regime studied in [Tl] corresponds to the case that each consumer node is surrounded 
by source nodes, since the installation cost is low compared with the transportation cost. 
When the installation cost is gradually raised, resource provision is achieved with less 
source nodes, but optimization requires the consumer nodes to be located in clusters 
surrounded by source nodes, forming the clusters observed in p^ . 

To formulate an algorithm, we introduce the constraints Si^i > for each deficient 
node. These contraints are not applied to surplus nodes as they are always satisfied. 
Introducing Lagrange multipliers Hi for the resource constraint, we minimize the 
Lagrangian 



with the Kiihn- Tucker conditions HiSi^i = and /Xj < 0. Optimizing L with respect 



particular set of {si}, we iterate these equations to find the corresponding set of {/ij}. 
The set of optimal {sj} is found by an approach similar to the the GSAT algorithm [31] . 
by comparing the Lagrangian in Eq. (j3]) for each choice of {sj}. In each step of this 
algorithm, a cluster of NfHp nodes is randomly selected. The network energies involving 




(2) 




(3) 





Figure 2. Simulation results of average energy per node and the fraction of network 
nodes acting as source nodes. Parameters: K = 3, (f>d = 0.9, N = 100, A'flip — 4, 100 
samples and 1000 flips. New clusters formed on increasing are sketched at the top, 
with filled and unfilled circles representing consumer and source nodes respectively. 



the different configurations of this cluster are compared, and the cluster configuration 
is updated to the one that yields the lowest network energy [H]. 

2.2. Major Simulation Results 

As shown in Fig. |2] for K = 3, two phases can be identified: (1) all-source phase for 
u^^ > a/3, in which all nodes are assigned to be source nodes due to the very high 
transportation cost; (2) partial- source phase for < u^^ < y/S, in which only some 
nodes are assigned to be source nodes. (In [13] we also identified a phase transition at 
= to an all-consumer phase.) 

The fraction of source nodes is a discontinuous function of u^^, showing abrupt 
jumps at threshold values of u~^. The step size of the curve decreases as increases, 
and gradually becomes unresolvable by the numerical experiments. This resembles the 
Devil's staircase observed in the circle map and other dynamical systems [26]. These 
threshold values of mark the positions at which certain configurations of the source 
and consumer nodes become energetically stable. Similar features are observed in RFIM 
due to the formation of ferromagnetic clusters resultant from the competition between 
the strengths of couplings and random fields [201 121] . Except for a shift of the average 
energy per node, these features are qualitatively similar to the simulation results of [H] 
which correspond to the case (/(^ = 1, if the latter is reinterpreted from the perspective 
of the source location problem. 

Measuring the average maximum cluster size of the consumer nodes in the samples, 
we observe abrupt jumps of the cluster size at the same threshold values. This indicates 
that new types of clusters are formed at each jump, as sketched in the top of Fig. |2j 
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Figure 3. Cluster of (a) source node, (b) singly consuming node, (c) doubly 
consuming nodes, and (d) triply consuming nodes. 



The observed threshold values can be calculated by considering the energies of consumer 
clusters surrounded by source nodes as shown in Fig. [3l obtained by the minimization 
of Eq. ([T]). By comparing the energy of different configurations, we have 

So > Si (singlet), 

So + Si>S2 (doublet), (4) 
So + 2Si>S3 (triplet), 

resulting in the threshold values in Fig. [2l These results agree with those obtained 
through the cavity approach in Section |3l We call the regime a/3/2 < < \/3 
with isolated consumer nodes the singlet regime, and v^21/25 < u'^ < the 
doublet regime. The isolated nodes are referred to as singly consuming, while the paired 
consumer nodes doubly consuming. 



3. The Replica Symmetric Ansatz 

3.1. The RS Recursion at Zero-temperature Limit 

We apply the cavity method [21 [3] assuming that the network has a locally tree-like 
structure. We denote as Ej{yj) the energy of the tree terminated at node j in the 
absence of its ancestor node i, when a current yj is drawn from j to its ancestor. 
Relabeling the descendents of j as k = 1, . . . , K — 1, Ej{yj) is expressed as 

E,{yj) = n{Ek=u Ek-i; Aj, yj). (5) 

The functional "H is given by 

n{Ei, . . .,EK-i;Aj,yj) 

2 2 

^min[ 5] E,{y,) + ^Q[-A,- + 2/.) + |] • (6) 

keCj\{i} keCj\{i} 

In the absence of node i, there is no supply or demand of resources through the cavity 
and the last term ?/j/2 should be absent. However, the presence of the extra term results 
in a clear interpretation of Ej, as we will see in the following sections. Care has to be 
taken when dealing with the change of the cavity energy, where yj is taken to be zero 
to eliminate the effect of the extra transportation cost on the dangling bond. 
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We note that Ej is an extensive quantity that depends on size of the tree. To 
formulate a recursion of an intensive energy, we write Ej{yj) as a sum of two terms, 

E,iy,) = Ejiy,) + E,iO). (7) 

We call Ej{yj) the cavity energy functions which correspond to the cavity fields in the 
language of the cavity approach, and represent the energy variation from Ej{0), as yj 
varies. In this case, Ej{0) = 0. Ej{0) corresponds to the energy of the tree when no 
current is drawn from the vertex. We further define the energy change AEj due to the 
addition of a vertex, 

AE, = E,{0)- J2 ^fcW (8) 

keCj\{i} 

which simplifies Eq. ([5]) to 

Ejiy,) = H«=i, • • • , E^-i, A„ - AE, (9) 

where 

AE,='H{EY,...,E]^^,;A„0). (10) 

We have thus separated the energy contribution due to the addition of a new vertex 
from the energy variation due to the changes in the current drawn from the tree. 

The distribution V[E^] of E^ over the vertices of the tree is given by the solution 

of 

K-l 

V[Ej] = / dA, n / dElV[E^]5[Ej - • • • , E^-i, K 

k=i 

+ ae{eX,...,eI_,-\,)]. (11) 

To elucidate the physical behavior of the system, we consider a node fed by K trees 
forming a Bethe lattice. For instance, we consider the average energy per node. The 
change in energy due to the additional node i is given by 

A^node = H(E]",...,E^;A„0). (12) 

Similarly, we can consider a link bridging two trees forming a Bethe lattice. The energy 
change due to the addition of a link between nodes L and R is given by 

A^Hnk(i?r, El) = min \El{y) + E^^-y) - ^1 . (13) 
y I Z 1 

Denoting (...) as the average over the capacities, the average energy per node is given 

by 

(A^) = (A£:,ode)-f (A^nk). (14) 
3.2. The Piecewise Quadratic Solution 

Due to the quadratic form of the transportation cost assumed in Eq. ([1]), we propose 
that the cavity energy functions are continuous and piecewise quadratic, namely, 

E^{y,)=mm[ft^{y,)]. (15) 
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Figure 4. An example of E^{y) composed of three quadratic fmictions (y) labelled 
by nk — 0, 1, 2. Each composite fmiction is characterized by its minimum position y^^, 
minimum value dj^ , and curvature . 



where rik = 0,1,2,.... Indeed, the recursive nature of the quadratic cavity energy 
functions have been fully employed in deriving the message-passing approach in [Qj [10] . 
We note in passing that a similar recursive structure was used in the Gaussian Belief 
Propagation algorithm [27] and apphed to processing continuous signals such as those 
in CDMA multiuser detection p8j. As a step forward, the ansatz in Eq. ( fTSl) further 
captures the multi- valley features in E^^y^), which is crucial in formulating the cavity 
messages for the present model. 

We call /^^ the ra^-th composite function of Ekiyk)- For Uk > 
quadratic function of the form 

whereas for = 0, to takes the form 
foiVk 



0, 



fLivk 



IS a 



(16) 



2 

y + Cfe + ye(yfe - ak) 



(17) 



The form of /o is relevant when node /c is a source node. Though /q is discontinuous, 
we will show that the resulting E^{yk) is continuous since the discontinuity at yk = cxk 
is masked by other quadratic functions. An example of a cavity energy function E^{yk) 
composed of three composite functions is shown in Fig. HI 

We denote the /^^ with the minimum d!^^ among all composite functions in E]^ as 



arg ( min d. 



is particularly relevant when we evaluate the energy of the system in the case that 
the ancestor node of is a source node. In this case, the resource of the ancestor can 
be freely drawn by node k without any consequences to other parts of the network, and 
the optimal current yk takes the value y^, . 

As will be shown in the next subsection, only a few composite functions are relevant 
in the singlet and the doublet regimes. Each composite function has its fixed values of y^^ 
and a^^ independent of k, but the constant terms is fc-dependent. Thus, the functional 
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n 


Vn 










i 

— 1 — 


1 




2iK-l) 


2 




- 1 


2{K^-K -I) 



Table 2. The table of cavity composite functions fn{v) = {VmO-n) from n = to 2. 



form of E^{yk) is effectively parametrized by tlie constant terms of the composite 
functions, as given by 

E^{y,) = {4,dldi...). (19) 

The relevant composite functions in the singlet and doublet regimes are shown in Table 

m 

3.3. The Recursion of 

Using the piecewise quadratic ansatz in Eq. f |T5|) . the recursion of E'^ in Eq. (jS]) becomes 
Ej{y,) = minH(/„t\ • • • , fn'^.^^, - mm AE,(/„\=\ . . . , f^~_\;A,). 

{nk} {nk} 

(20) 

We first consider the recursions for a deficient node j. Suppose the node is assigned 
to be a source node, i.e., + J2k€Cj\{i} ~ Vj < 0- this case, 

= f + y®(^^-^^- C)+ E <' (21) 

keCj\{i} keCj\{i} 

which is exactly the form of /o in Eq. f lT7|) . The combination {n^,} of composite function 
minimizes the last term in Eq. f l2T|) . Subject to a vertical shift by AE'j, Eq. fl2T|) is taken 
to be /q in Eq. f[T7|) characterized by 

^^o = y+ E ^n*-A^. (22) 

fcg£j\{i} 

and 

fcez;j\{i} 

For the moment, we ignore the possibility that combinations other than {n^} may result 
in further discontinuity in /q, as we will show that composite functions /„ with > 
eliminate the effect of the discontinuities. 
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Next we suppose the node remains as a consumer node. In this case, Aj + 
T.kec,\{i}yk -Vj = 0' and mm{n,}U{f^=\ . . . J^-_\; Aj^yj) is computed subject to 
this equahty constraint, 



keCj\{i}yk-yj=^ 

k&Cj\{i} 



(24) 



where 



1 + 



1 



Y, 



{nk} 



1 + E 



keCj\{i} 



(25) 



D 



{nk} 



1 + E 



keCj\{i}i'^^nk 



and the optimal currents drawn from the descendents are 



Vk^iUk} 



«e£j\{j}(20'nj 



+ yk- 



(26) 



Next, we consider the recursions for a source node j. In this case, the cavity energy 
function consists of the composite function with Uj = only, with 

di= E ^n;-Ai?,. (27) 

Table [3] summarizes the combinations {n^} which lead to the composite functions 
= 1, 2 in Table [21 enabling us to analyze the singlet and doublet regimes. We represent 
these relations of {uk} with nj by the mapping Ai. As an illustration, the relations in 
Table E] can be expressed as 



A^(0,...,0) 



M(1,0, 



,0) 



f28) 



The constant term dl is taken to be the minimum constant term in all combination of 
/^^ with M{{nk}) = Uj, yielding 

di. 



mm 

{{nk}\nj=M{{n^.}) 



D 



dt 



-AE, 



(29) 



Remarkably, the functional recursion of in Eq. (|9]) is now simplified to a recursion 
of the constant terms in Eqs. (122]) and f l29p . 

The physical interpretation of the composite functions Uj = 0,1,2 is revealed by 
considering the patterns of optimal currents around node j. For nj = 0, no resources 
are drawn from the ancestor, i.e., yo = 0. Moreover, the optimal currents drawn from 
the descendents are fjn*, which are non-positive according to Table |3l This shows that 



Uj = corresponds to a resource providing state (see Fig. El^a)). 
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{nk} = 
(ni,n2, . . . ,nK-i) 


Y{rn.} — Vrij 










(0,0,. ..,0) 


1 

~K 


K 

2{K-l) 


1 

2K 


(-- -) 


1 


(1,0,. ..,0) 


1 


- 1 


(K+l) 




2 


K -1 


2{K^ -K-l) 


2K{K-l) 


(°'ir-l''"'ir-l) 



Table 3. Combinations {nk] from descendents that lead to the composite functions 
rij = 1, 2 given by Tabled Permutation of (rii, 712, ... , tt-a'-i) result in the same rij. 




Figure 5. An example of minimization by 'H{f^~-^,...,f^, ^■,Aj,yj) with the 
combination {n^} = (n^, 71,2, . . . , The labels Over Satisfied, Just Satisfied, 

and Unsatisfied refer to the conditions > 0, = and < respectively. The 
just-satisfied portion of H intersects with /q at yj — aj as shown. 



For Uj = 1, currents of are drawn from the ancestor and all descendents, i.e., 
yl^i = . . . = y^-i = K~^, as shown in Table [3l Hence nj = 1 corresponds to the singly 
consuming state, (see Fig. [3]^b)). 

For Uj = 2, currents of {K—l)~^ are drawn from the ancestor and K—2 descendents, 
leaving the link to the remaining descendent idle. This corresponds to the doubly 
consuming state (see Fig. EJ^c)). Nodes with uj > 1 correspond to other modes of 
resource consumption. 

We return to verify that is continuous piecewise quadratic, despite the 
discontinuity of /q in Eq. ( fTTl) . This can be done by substituting yj = aj in Eq. (l23ll 
into the composite function Eq. (IMl) characterized by ^{n*}, ^{n^} ^iK} ^'i' ^ 1^ - 
The result shows that T-L = a|/2 + J2k€C ^n*- Hence it intersects /J at the lower end 
of the discontinuity. As shown in Fig. |5l the discontinuity is masked. 

We finally derive a simplified expression of AE'j from Eq. ( [T0|) . Since 
^(/ni' • • • ' fuK-i'^^J^ 0) is already obtained in Eqs. (ED and (El]) with yj = 0, 



AEj = min 



keCj\{i} 
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^0(-A,)+ J2 <]• (30) 

keCj\{i} 

4. Closed Sets of Cavity Energy Functions: The Intense Simplifications 

4.1. The Intense Simplification 

In this subsection, we consider networks with no surplus nodes (0s = 0). In the singlet 
regime, the energetically stable configurations consist of only the source nodes and singly 
consuming nodes. Hence we consider ii^^with only /q and /i as composite functions, 
is thus given by 

EUyk) = {4,d',), (31) 

as a simplification of Eq. ( [T9|) . Composite functions with > 2 have > d^* in this 
regime, and their corresponding configurations are not stable. The recursion relations 
of do and di in Eqs. (122]) and ( 12^ are simphfied to 

di = Y+ E min(4dt)-AE,, (32) 

keCj\{i} 

d{ = 7^+ E d'o-AE,, (33) 



2K 

keCj\{i} 



where 



AEj = min 



^ min(4,c/t),^^^— ^+ E rf^ 



(34) 



k€Cj\{i} ^ ' fcG£j\{j} 

To determine the pattern of current fiow, it is sufficient to consider the recursion of 
ej = d{ — d^, given by 

^j = -l- Yl "^in(0,efc), (35) 

keCj\{i} 

where 

1 

7 = . 36 

The simple recursion leads to a closed set of K cavity energy functions 

E^{y) = {di,di + q^) withg = -l,0,...,ir-2. (37) 

These functions are classified to be consuming for q = —1, bistable for q = 0, and resource 
providing for q = 1, . . . , K — 2. Their absolute minima are located at y = —K^^, both 
y = and —K^^, and y = respectively. We call their states c state, b state and s state 
respectively; examples for the case of = 3 are shown in Fig. |6l The b state behaves in 
the same way as the s state in the recursion relation, but physically they correspond to 
different cavity states. The differentiation between b and s states is required only when 
the entropy of the ground state is calculated. For most other purposes, grouping b and 
s states together further simplifies the analyses. We thus denote the b and s states as 5" 
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Figure 6. The closed set of Ek{y) eA K = and u ^ 5/3, corresponding to (a) the 
c state , (b) the h state and (c) the s state E^ ■ 
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Table 4. The recursion relations given by Eq. The column AEj is valid for 

u^^ > \J K ~ \ in the singlet regime; a constant {K — l)[u^/2 — \j2(K — 1)] has to be 
added to Ai?j for < \JK — \ in the singlet regime. 



state , and the c state as C state in subsequent analyses. Their recursion relations are 
summarized in Table HI (The source state (s state) in this paper should not be confused 
with the satisfied state (also denoted as s state) in [II]. In fact, the s, h states in [H] 
have the same cavity energy functions as the c, s, h states in this paper respectively.) 
The recursion relation in Table H] can be summarized by the symbolic equations 

^+... + S^C. (38a) 

K-l 

all other combinations — )■ S. (386) 

To calculate the average energy per node, we can apply the same simplification to A£^node 
and Af^iink in Eqs. f|T2|) and f lT3|) respectively. We denote the jull states of a node as 
C, B and iS, respectively representing the consuming, bistable and resource providing 
states. They are obtained symbolically via 

S+... + 5^C (39a) 

K 

P + g + . . . + g ^ B (395) 

K 

all other combinations — )■ 5. (39 c) 
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The energy changes are 



A^nodc = y + min 



K 



(40) 



A£:iink = min(0,eL,eR). (41) 

These expressions are vahd for > y/K — 1 in the singlet regime; Nc[u'^ /2—1 /2{K —1)] 
has to be added to the expressions for u^^ < \/K — 1 in the singlet regime, where Nc is 
the number of vertices in the C state . 

Numerical iterations of Eq. f l35|) starting from random ek show that the closed set 
of is stable. The closed set corresponds to the integer cavity fields in the language 
of the cavity approach. Another example of a closed set of is found in the doublet 
regime as described in Appendix. 

We return to discuss the range of the singlet regime. From Eq. (1351) . it becomes 
apparent that when 7 < 0, we would have ej always positive, implying that the singly 
consuming state is always unstable. Hence a necessary condition of the singlet regime 
is 7 > 0, or < \/K. 

Similarly, in the doublet regime, we obtain the recursion relations 

d{-di= --f- J2 min(0,rft-4,rft-4), (42) 

k&Cj\{i} 

dl- dl= - K + mm{d'l - d^) - min(0, d'l - d^, d'l - d^), (43) 

fce£j\{i} 

where 

Consider the difference d2 — d{ = —k + min{di — d^) + 7. In the singlet regime 
described by Eq. fl35l) . d^ — d\ is always positive only if k < 0, implying that 
u^^ > a/ K{K — 1)/ {K + 1). Combining the two results, the range of the singlet regime 
is ^^/I^(l^^^T}J(I{^^T) < u^^ < \/K, agreeing with the result reported in Fig. |2l 

It is convenient to represent the recursion relations in the probabilistic framework 
of BP algorithm [29j. We denote as the probability that node j is in the C state , 
in the absence of the ancestor node i. The probability that node j is in the S state in 
the absence of i is described as ipi^\ We call V'c^* cind ipi^^ the cavity probabilities. In 
this framework, Eq. (138 a\i in the singlet regime can be written as 

i^r= n (45) 

kec,\{t} 

and ipi^^ + V's"^* = 1- Note that the recursion relation (|5]) with the functional form 
of E^ is now represented by a recursion of a single scalar ip^c~^^, greatly simplifying 
the analysis. This framework will be useful in establishing a connection with the RSB 
ansatz. 
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4-2. Networks with Surplus Nodes 

In this subsection, we consider networks with surplus nodes (0^ > 0) in the singlet 
regime. The recursion relations of the deficient nodes follow those in Table IH whereas 
the of the surplus nodes are always in the s state. The symbolic equations Eq. (138 al) 
are thus extended to 

+ S^^^C. (46 a) 




all other combinations ^— 7- ^ S*. (46 &) 

all combinations S. (46 c) 

The recursion of ipi'^^ is given by 

^r = 5A„-i n (47) 

keCj\{i} 

and tpi-^' + = 1. 

The full energy change is given by 



■^^node 



+ mm 
2 



A,- = -1 



(48) 



^min(0,e{), Ai = A 

and Af^iink is still given by Eq. ( I40l) . For deficient nodes, the optimized state is either 
consuming, bistable, or resource providing. Surplus nodes are fixed to be source nodes. 
A^node and the full states of node i are shown in Table [5] with the corresponding 
combination of descendent states. The full states of the node i are thus described 
by the probabilities ipc^ i^h given by 

^B = 5A.,-iJ2^t' n (l-^c^O, (49) 
je£i k&Cj\{i} 

with ip^, + ip^Q + ipg = 1. 

We derive the fraction of nodes with different full states by assuming the 
independence between descendent branches of a tree. From Eq. ( H5l) . the average 
probability (ipc) of a node in the cavity C state is thus given by 

{^p,) = M^-{^c)f-\ (50) 

where (. . .) represents averaging over nodes. Iteration of Eq. (H5l) on Cayley trees 
reveals that (tpc) does not approach the stable fixed point of Eq. fl50|) when (p^ is high. 
At = 1, a period- two solution of (tpc) = and 1 emerges. Physically, this corresponds 
to alternating layers of consumer and source nodes on Caylay trees. This happens in 
similar problems such as the Bethe glass |23j and is referred to as the modulation mode. 
On real networks, nodes are randomly connected, rendering whole layers of consumer 
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and source nodes highly unlikely. The period-two situations are suppressed, making a 
fixed point solution of Eq. f l50|) possible in random networks. The fraction of nodes with 
full states C, B and S are thus given by 

/c=0,(l-(^.))^, 

fs = Kct),{^,){l-{ij,)f-\ (51) 
with fc + fB + fs = 1- We leave the discussion of their physical interpretation to Section 



5. The Average Energy 



We evaluate the average energy by considering the energy change due to the addition 
of new nodes and links. Summarizing Table [5], (A£^node) becomes 

;i-(V^,))^7l -^(V^c)7. (52) 



u 



(A^node) = (Pd 

Similarly, {ASunk) is given in Table |6] by 

After some algebra, Eqs. (JH]), ([50]) ([52]) and ([53]) lead to 



^1 



u 



RS 



7- 



(53) 



(54) 



When a deficient node changes from a source to a consumer, the energy of the node 
reduces by 7 from We thus identify the coefficient of 7 in Eq. flM)) to be the 




Figure 7. The average fraction of source nodes in simulations with K — 3, 
N = 100, 500 and 3000A^ steps with iVflip = 4, as compared with the ff^^ obtained 
from the RS ansatz. Inset: ff ^ for different values of K. The circles indicate the 
values of 0^ above which the RS assumption is not stable. 



fraction of consumer nodes. The fraction ff- of source nodes is then given by 

/f = 1 - i^c) - (f - l) (^^)'- (55) 

Note that ff-^ is distinguished from ff^, since fs = fs + /b/2, i.e., fs also counts those 
bistable nodes that become source nodes in the network configuration. Through the 
linear relationship {S) = (j)dU^/2 — /^^7 = 1/2K — (1 — 0d)u^/2 + we consider fg 

as a measure of the average energy (S). 

The inset of Fig. [7] shows the fraction of source nodes as a function of (pd derived 
from Eqs. fl50|) and fl55l) . For all connectivities K, ff-^ decreases with (pd- A higher 
connectivity leads to an increase in ff'^ since more nodes are required to convert to 
source nodes to satisfy the demand of a consumer node. 

Figure [7] shows the difference between ff-^ and /|™ obtained from numerical 
simulations, in which the energy of real instances is minimized by the GSAT algorithm 
as described in Section [2l2l The differences between ff'^ and /f™ are roughly zero when 
(j)d is below some critical value. Above the critical value, ff'^ is significantly lower than 
/I™. Hence the energy is lower than the simulated energy. This discrepancy is 
related to the instability of the RS ansatz, which will be discussed below. 



5.1. The Soft Nodes in the Ground States 

Nodes with the full state B are bistable between the consumer and source states. 
Figure |8] shows an example with two central nodes in the i3-state. Using Table HI 
the cavity energy functions passed among them are worked out. In particular, is 
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(a) 




Figure 8. (a) An example of a configuration with two central soft nodes, annotated 
with cavity energy functions passed among them, (b-c) The two degenerate states of 
configuration in (a). 



sent from the bistable nodes to all its resource-providing neighbors, implying a zero 
current or a current of —K~^ would have no effect on the optimized energy. Hence, the 
central pairs correspond to the soft nodes which can be a consumer or a source node in 
different degenerate configurations. These degenerate configurations are connected in 
the configuration space, since transitions among them only involve the flipping of states 
of the soft nodes. The RS ansatz assumes that the configuration space is dominated by 
a single cluster whose states can be obtained from each other by flipping the states of 
the soft nodes [30]. In other words, using the notions of S, C and B states, the single 
cluster is described by a unique set of iS, C and B labels for each node. On the other 
hand, if the network enters the RSB phase, the configuration space is dominated by 
numerous clusters, each with its own set of 5, C and B states of the nodes. Transitions 
among the clusters involve flipping the hard (non-soft) nodes as well. In this section, 
we discuss the RS case. 

The fraction of soft nodes is given by the fraction of nodes with full state B. From 
Eq. (ISTl) . we have 

ffl = fB = KUi^^){l-{i;^)r-\ (56) 

On the other hand, the hard nodes are either consumer or source nodes in all degenerate 
states. They correspond to nodes in the backbone in a vertex cover [H [71 [8]. Surplus 
nodes are certainly in the source backbone. Deficient nodes can be found either in the 
consumer or source backbone. From Table [5l the fraction of nodes in the consumer 
backbone is given in Eq. flSTl) by the probability of finding a node with full state C, 
namely, 

/con-bone = /c=0d(l-(^c))^ (57) 

Remarkably, it can be shown easily /so|/2 -|- /^n_bone = ^ ~ ff"^ which implies that 
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Figure 9. The average fraction of soft nodes in simulations of A' = 3 and 3000 
steps with iVfiip = 4, as compared with f^%- Inset: ff'^/(f>d for different vakies of K. 
The circles indicates the values of 0^ above which the RS assumption is not stable. 



exactly half of the soft nodes are consumers. 

The fraction of soft nodes relative to the fraction of deficient nodes obtained from 
Eq. (l56l) is shown in the inset of Fig. |9l When the fraction 0^ of deficient nodes increases 
from 0, the fraction increases until it reaches a maximum value at (ipc) = and 
(f)d = K^~'^/[K — When (pd is small, most deficient nodes are surrounded by 

surplus nodes and hence are found in the consumer backbone, leading to a small fraction 
of soft nodes. When 0^ increases, the probability of finding contiguous deficient nodes 
increases, which leads to an increase in 

The fraction of soft nodes in simulations is compared with in Fig. [91 In 
simulations, we use an algorithm similar to the GSAT [31] to lower the energy of the 
system until it becomes steady. Then we allow further cluster flips and identify the 
nodes which can be flipped with no change in global energy, such as the example in Fig. 
m^b) and (c). Similar to Fig. [TJ the simulation results have an excellent agreement with 
the RS result when 0^ is below the value at the peak of f^ft/^d- Above this critical 
value, the fraction of soft nodes in simulations is less than the RS prediction, which will 
be shown to be due to RS instability. 

Algorithmically, a high fraction of soft nodes leads to an unfavorable consequence. 
After the convergence of BP on real instances, nodes with B state are bistable and 
a further determination of the final optimal configuration is required on the resulting 
sub-graph of B nodes. A random assignment of B nodes to either consuming or resource- 
providing does not generally result in an optimal configuration. When (p^ is low, the 
sub-graph of B nodes are disconnected and the assignment is easy. When 0^ is high, 
the sub-graphs of B nodes are connected and the assignment is more difficult. 
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5.2. The Instability of the RS Ansatz 



The discrepancy between the simulation results and the predicted average energy S-rs 
and /g^l suggests that the RS ansatz is unstable at high 0^. In the RS formalism, a 
single ground state is assumed. We thus examine the stability of this assumption against 
the picture of multiple ground states, by relaxing the constraints of = 0,1. In other 
words, < ■ipi'^^ < 1 
variables 



To study the stability of the integer ansatz of V'c^*; we define the 



and r/^^* 

integer values of ip; 



1 _ ^.-^ _ 



I a 



v'r = 5^rM (58) 

Thus, the recursion rule f l45|) is extended to include non- 
1 corresponds to the onset of non-zero probabilities 
over the range of < ip^^"^ < 1, indicating the occurrence of glassy behavior with 
node j having a probabilistic distribution of C and S states in the stable states of the 
configuration space. Hence the RS solution is stable only if r/^"^* = for all j — )■ z is a 
stable fixed point. 

We thus formulate the recursion relations of ?7^^*, rji^^ and r]^"^* by considering 
Table HI namely, 

K-l 



1 n 



111^, 

k=l 

K-l 



K-l 



k=l 



IK". 



k=l 



K-l 
k=l 



k^j 



(59) 



as (rjc), {rjs) and (rjg) respectively, we 



Denoting the site average of r/^"^*, 'ql^^ and 
obtain their stable fixed points by solving 

ivc) = Mvsf-' 

ivs) = I - Mivs) + {%) f-' 

(Va) = MiiVs) + {r,,)f-' - {r,;)^-'] (60) 

For all (j)d, {rig) = is a trivial solution of the last line of Eq. (l60ll . and the above 
recursions reduce to the RS recursions fH5|) . By introducing a small perturbation Srj^^^ 
to rjg^^ = 0, we obtain the corresponding Sri^g~^\ The solution of (rjg) = is stable under 
the perturbation if 



{K-l)MVs) 



K-2 



< 1, 



(61) 



Alternatively, the stability of the RS solution can be studied by considering the 
propagation of fluctuations ((^V'c"^"')^) under the recursion relation Eq. (H5ll [31] . This 
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Figure 10. (a) The stable solution of Eq. (|60)) for K ~ 3. The symbols represent the 
simulated fraction of non-converging BP messages, (b) The dependence of the variance 
on (pd- The dotted lines indicate the corresponding values of (p^^ ■ 



leads to the Almeida-Thouless (AT) stability condition, 



2\K-2 



< 1. 



(62) 



In the RS regime, ((V's)^) = (Vs) since ipl^"^ = or 1. The AT stability condition is thus 
equivalent to Eq. ( I6T!) . 

The stable solution of (r/c), (Vs) and {rjg) foTK = 3 is shown in Fig. [TOTa). The RS 
solution becomes unstable when 0^ > 0.75 for K = 3. This critical value agrees with 
those found in the simulation of the average energy and the fraction of soft nodes. For 
general values of K, simple algebra leads to the AT line 

which separates the RS and the RSB phases in the space as shown in the inset of Fig. [TTl 
From Eqs. (ED]), and (ESD, it can be shown that {tp^) = /sol = on the AT line. 

In the large K limit, 0^""" approaches e/K. This result has an interesting connection 
with the vertex cover problem. Considering the covered set as the set of source nodes, 
all links involving surplus nodes are covered. The remaining links are those among 
the deficient nodes. These deficient nodes have at least one neighbor being a deficient 
node. Hence from Table [5l their surplus node neighbors do not affect the states of 
the deficient nodes. Rather, their states are determined by the states of their deficient 
node neighbors. Thus, the problem of minimizing the covered set size reduces to one 
that minimizes the subset size of covered nodes in the subnetwork of deficient nodes as 
sketched in Fig. [TTl In the large K limit, this subnetwork has a Poissonian connectivity 
distribution with a mean K(j)d. The result Kcj)^^ = e agrees with the point of RS 
instability derived in [6l[7]. 
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Figure 11. Sketch of the Poissonian subnetwork of deficient nodes in the large 
connectivity limit with only a fraction 0{K~^) of deficient nodes (squares) and the 
rest being surplus nodes (circles). Filled and unfilled symbols represent consumer or 
source states respectively. Inset: The RS and RSB phases in the K-(j>d space. 

Comparisons between the simulation results and the RS analytical results from Figs. 
[7] and [9] have yielded evidence of an AT transition in real instances at ipd = ■ Here we 
provide two more supporting numerical experiments. In the first experiment, we consider 
the BP algorithm initialized with ipi~''^ = 0, 1 for all j — )■ i. As shown in Fig. [TOT a). 
effectively all messages converge to their steady states in the RS regime. However, a 
significant fraction of messages fluctuates between and 1 when 0^ rises above cj)^ , 
indicating the breakdown of the RS ansatz. It leads to the non-convergence of the BP 
algorithm on real instances. As shown in Fig. [roT a). the fraction of non-converging BP 
messages has an excellent agreement with the theoretical values of {rig). Consequently, 
improved algorithms are needed for ground state searching. Algorithmically, decimation 
procedures, such as those used in the survey propagation (SP) algorithm [4], are required. 
We have tested the BP algorithm with decimation and obtained results with simulated 
energy lower than the GSAT algorithm, which will be reported elsewhere. 

In the second experiment, we consider numerical iterations of Eq. f H5l) using 
population dynamics [9l[T0]. We start with different sets of initial values of = 0, 1, 
each set following the same sequence of random connections among the nodes. When the 

dynamics reaches the steady state, we measure the variance Cj^i = (V'c"^*)^ "~ (^'^c~^*j > 
where the overline denotes the average over random initial conditions, crj^i ~ implies 
that the cavity state of node j is independent of the boundary condition of the tree 
represented by the population dynamics, and is always frozen in either S* or C states. 
(Tj^j > implies that the cavity state of node j is dependent on the boundary conditions 
and shows a long range correlation. The numerical results of o'j-^i averaged over nodes 
are shown in Fig. [TOT b) for K = 3,4,5. In the figure {o'j^i) > when (p^ > (p^^ with 
almost no dependence on the number of iterations in the population dynamics. 
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6. The One-step Replica Symmetry Breaking Ansatz 
6.1. The IRSB Formalism 

In the RSB formalism, the network behavior is explained by the dominance of numerous 
states in the configuration space, instead of a single ground state in the RS formalism. 
Here, we consider the IRSB ansatz [321 [33] where the density of states with energy 
e = S/N per node is assumed to have the form A/Ar(e) = exp[A^E(e)], for a network of 
nodes and total energy S. S(e) is called the complexity or configurational entropy. 
For small changes in the average energy, we can write S(e) = a;(e — e^^) where e^^ is 
the reference energy. This assumption leads to a recursion for the cavity probability 
functional Pj [E^] 



PAEj{y,)] = ^ n {j'^El{y,)P,[El{y,)] 



keCj\{i} 

X W {5[Ej{y,) - U{E\, . . . , El_,- A„ y,) + /\E,{E\ , . . . , El_,- A,)] 
yj 

X exp[-xAE,«, . . . , El_,- A,)]) (64) 
where S is the normalization constant. We now define the right hand side of Eq. (16^ to 



X 



be Hp and the recursion can be represented by Pj[Ej] = 'Hp{Pk=i, . . . , Pc-i, Ej; A 
Solving the recursion in Eq. f l64|) by population dynamics is equivalent to solving for a 
stable functional distribution Q in 

K-l 

Q[P,] = / rfA,p(A,) n / -DPkQlPk] (65) 
k=i 

X H 6{P,[Ej] - 'Hp{Pk=i, Pc-i, Ej; A,, x)}, 

EJ 

which is an analogy to the RS case of solving for in Eq. ([IT]). Note that the 

RS recursions in Eq. ^ correspond to the recursions of the cavity energy functions 
E^ yielding a solution of the functional P[E^], while the IRSB recursions in Eq. (164)1 
correspond to recursions of the functional probability P[i?^] yielding a solution of the 
probability functional Q[P]. 

To analyze the physical properties of the network, we write the partition function 
when K cavity probability functional feeds a central node. 

Hnode = n ( / T^^lPAE^]) exp[-xA W(^r, ...,El- A,)]. (66) 

with A£^node giveu by Eq. f [T2]) . Similarly, the partition function obtained by bridging 
two trees with a link is 



--link 



j VEIVEIPl\EI]Pr[EI] exp[-xA£H„k(i^;r, E^)]. (67) 
The average configuration free energy is given by 

^x) = -^((InSnodc) - y(lnSii,k)) (68) 
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The averages (. . .) are taken over Aj and Pj from the distribution Q[-Pj]. $ is related to 
the complexity S and the average energy e by 

x^ = xe- S, (69) 

with dTj/de = x. S and e are parametrically dependent on x via 

E = x^^, (70) 



dx 
d{x^) 



(71) 



dx 

While solving for P[E^] in Eq. (|64]) is in general difficult, simple solutions can be 
obtained if a closed set of countably many is sufficient to describe the recursions of 
. The singlet regime in which the S and R states form a closed set is a good example. 
We emphasize, however, that the techniques are generally applicable to regions beyond 
the singlet regime where closed sets of E^ are found, such as the commensurate point 



in the doublet regime discussed in Appendix A 



6.2. The IRSB Solution 

In the singlet regime, there are only two representative states, C and S states, we 
parametrize P[ii^^] as 

PAEJ] = - E^) + mEj - E^). (72) 

Using Table H] to obtain /S.E for different combinations of C and S states, Eq. ( l64l) can 
be simplified to 



k&C,\{i} 



+ 5/ 



^r=5A„-i^7^[ n (^c^^e^" + ^r^') - n 

fcG£j\{i} fcG£j\{j} 

The above recursions of V'c"^* are gross simplifications of the recursions of the functional 
probabilities in Eq. fl64|) . These equations can be solved by population dynamics 
involving a pool of values of '?/'^^' with < V^^^* < 1. Alternatively, the solution 
to the IRSB recursion can be found by directly solving for the distribution Q{ipc) in 
Eq. (j65l) which is isomorphic to Q[P] in Eq. (j65l) . 

The stable solution of Q{ipc) for K = 3 and x = is shown in Fig. IT^ a). When 
4>d < 4>d^y there are no fractional components of ipc, and Q{ipc) = {4'c)S{4'c ~ 1) + (1 ~ 
{ipc))S{ipc)- When 0^ > (j)^^, non-zero components of < ipc < ^ exist. This agrees with 
the result in Section I5l2] that in this regime, the RS solution is an unstable solution of 
Eq. f l65p . When 0^ = 1, there is no disorder in the capacities. All vertices are identical 
and Q{ipc) = 5(^c — (V'c)) where {ipc) is given by the RS equation fl50|) . This means that 
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Figure 12. (a) The stable solution of Q{ipc) obtained by solving Eq. (j65|) in the 
singlet regime using Eq. d75|) with K — 3 and x — 0. (b) The stable solution of Q{ipc) 
with K = 3 and (pd — ^■ 



among the different states of the system, all vertices are equally probable to be in C 
state with probability (ipc)- The stable solutions of Q{ipc) are dependent on x via the 
factor x'-y. In the case of 0^ = 1, the dependence of Q{ipc) on x is shown in Fig. [T^b). 
The position of the delta peak at {ip^) = (1 - (V^c))^"V[(e'^^ - + 1]^"^ shifts to 

left from X7 = to 0:7 = x*^^'*7 = 1.92. For x > x*-^-* in Fig. [TW b). Q{ipc) becomes a 
continuous distribution, indicating the instability of the IRSB ansatz to be discussed in 
the next subsection. 

With the stable solution of Qii'c), we derive the complexity E(e) in the singlet 
regime. Obtaining A^node in Table |5l we write the partition function Snode in Eq- dM]) 
as 

-xu'^/2 



• — ^node 



5^ 



-1 1 n 



c ^ 



,2:7 



1) n 



(74) 



With Ai?iink from Table |6l we write Snnk as 



--link 



3^7 



(75) 



The configurational free energy $ is given by Eq. ( 168|) with averaged over Qlipc)- 
The complexity S(e) obtained is shown in Fig. [131 Again, e is expressed in terms of the 
fraction of source nodes through the relation fs = [e — 1/2K + (1 — (pdju"^ /2]/'j. 

Generally, we identify three segments on the complexity curve: (i) the unphysical 
segment (the dotted segment with S > 0), (ii) the physical segment (the solid segment 
with E > 0), and (iii) the highly unlikely segment (the dotted segment with S < 0). 
The physical segment of S is related to the number N' of metastable states with energy 
e. We denote the values of x when S = and S is maximum as Xg and Xd respectively. 
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Figure 13. The IRSB complexity I](e) obtained by numerically solving for Q{ipc) 
with K = 3 and 0^ = 0.94,0.96,0.98, 1. Inset: I](e) for K ^ 3 and (jja = 1 obtained 
from the ansatz Q{ipc) — S{ipc — (V'c))- It approaches the modulation limit at fs = 0.5 
as X ^ oo. Symbols are spaced at intervals of 0.9 in xj. 

e{xs) corresponds to the lowest energy among the states with non-vanishing complexity, 
which is considered as the ground state in the picture of IRSB. 6(0:^) corresponds to 
the energy of the states with the highest complexity, which is believed to be the states 
where search algorithms get trapped, giving rise to dynamical transitions. However, 
recent work on the coloring problem shows that the efficacy of the BP algorithm is not 
affected by the dynamical transition [30]. We leave this issue for future studies. 

The segment of negative complexity corresponds to states with vanishing number 
in the thermodynamic limit. Its physical meaning is clear in the limit x — )■ 00, which 
corresponds to a single state with lowest possible energy, since the reweighting process 
allows only one state. We show S in the inset of Fig. [13] as x — )■ 00. The result is 
obtained from the solution of Q{ipc) restricted to be the the IRSB solution in the form 
of Qii^c) = ^{i^c ~ {''Pc)), i-e. the unstable solution when x > x^^\ The complexity curve 
approaches the limit of ipc = 0.5, which corresponds to the highly unlikely structure of 
networks with alternating layers of consumer and source nodes. 

6.3. Comparison with Real Instances 

Though successful in finding low lying states on real instances, the GSAT algorithm we 
described in Section 12.21 requires long computation time for large systems. To compare 
the predictions of RS and IRSB approximations with real instances, we employ the 
BP with decimation. As the BP does not converge in the RSB phase, we measure 
the time average of the fiuctuating messages V'c^* and evaluate for each node the time 
average probability of the full iS-states. The node with the highest ipg is fixed to be 
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Figure 14. The optimized fs on real instances as a function of A'^ with = 3, as 
obtained by the GSAT algorithm (o) and the BP algorithm with decimation (A). The 
horizontal dashed lines show the analytical results of the RS ansatz, fs{xs) and fs{xd) 
in the IRSB ansatz. Inset: the fraction /dec of decimated nodes before convergence. 



resource providing. Only source nodes are decimated, as bistable nodes should be left 
undecimated, and decimating a consumer node will fix its neighbors simultaneously, 
which may hinder the convergence of the BP. By repeating the above procedure, the 
BP messages finally converge and the full state of all nodes are determined. 

Figure [H] shows that lower simulated energy can be obtained from the BP with 
decimation as compared with the GSAT algorithm. The simulated energy approaches 
the e{xd) as obtained by the IRSB ansatz. However, the fraction of decimated nodes 
before convergence increases with system size as shown in the inset of Fig. [TH This 
greatly increases the computational time for large systems if only a single node is 
decimated at a time. Several nodes can be decimated simultaneously to shorten the 
computation time, with a tradeoff in energy. 



6.4- The instability of IRSB formalism 

To test the stability of the IRSB ansatz against further steps of RSB, we consider two 
kinds of instability that leads to the two-step RSB (2RSB) formalism [35]. They are 
the so-called type I and type II instabihties, corresponding to the aggregation of states 
and fragmentation of states respectively, as shown schematically in Fig. [151 In type 
I instability, metastable states aggregate in clusters and the IRSB ansatz is valid in 
each cluster. The whole state space is composed of clusters and the 2RSB formalism 
is required to describe the state space structure. In type II instability, some states 
split to form clusters of states instead of single states. The IRSB ansatz is valid inside 
the clusters while 2RSB ansatz is required to describe the state space structure. It is 
generally beheved that once IRSB is not stable, the full RSB is required to describe 
the system, as illustrated by the IRSB instability found in graph coloring [36] and K- 
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Figure 15. The schematic picture of the state space in the RS, IRSB and 2RSB 
phase. Black circles represent singly connected states. Grey circles represent clusters 
of states, k is an arbitrary vertex in the network. Its probabilities of the C states 
in the IRSB phase, and several clusters in the 2RSB phase, are shown in the figure to 
illustrate the aggregation of states. Its cavity energy E]^ in an arbitrary state in the 
IRSB phase, and several of its fragmented states in the IRSB phases, are also shown 
in the figure to illustrate the fragmentation of states. The values of ip'^ are assigned 
for the purpose of illustration only. 



satisfiability problems [37]. Here we focus on the case of i^' = 3 in the IRSB regime 
with (pd > (p^^, following the approach of |36] . 



6.4- 1- Type I instability: Aggregation of states In the IRSB phase, each vertex k 
is characterized by the cavity probability functional Pk[E^]. This characterization is 
much simplified due to the small closed set of states in the singlet regime, rendering it 
isomorphic to single values of ip'^. In the 2RSB formalism, each vertex k is characterized 
by the probability functional of P^- For the closed set of states in the singlet regime, this 
characterization is isomorphic to the distribution Qk{ipc)- example illustrated 

in Fig. [151 the probability ip'^ = 0.4 is originally uniform in the state space of the 
IRSB phase. When it enters the 2RSB phase, ip'^ starts to take up diversified values of 
ip^ = 0.39, 0.40, 0.41 in three different clusters of states. 

To examine type I instability, we test for the possible spreading in tp^"*^ by 
introducing small fluctuations in ip^^^ of the descendents. We define Tk^j{x) to be 
the derivative 

dtp^-"' 



dtp, 



IRSB 



(76) 



from Eq. (173|) . which is evaluated at the values of ip^^^ from the IRSB solution. The 
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propagation of noise is thus given by 

m-^r= E [Tk^A^Wis^'c^^?. (77) 

k&Cj\{i} 

In the thermodynamic hmit, we consider a tree structure with L generations. The noise 
at the L^^ generation vanishes if (Tl-i^lTl-2^l-i ■ ■ ■ Tq^i)^ — as L — > oo. Hence 
we measure the quantity 

Xi\x) = {K-lf{TL-l^LTL^2^L^l...n^l)^ (78) 

averaged over the quenched disorder and boundary conditions of the trees. This can 
be implemented by iterating Eq. fl73|) using population djTiamics and evaluating the 
corresponding X^l\x) in each iteration. Alternatively, we solve for the distribution 
VlI^^l ,'ipc] layer L by the recursion relation 



k=2 



K-1 



—xu /2 

0,5[aw , - (Ti^L+i)^A?]5[^^^ - ^^(1 - i^^) n (1 - 

k=2 



+ MXtUm:^^')\ (79) 
with the initial condition 

P,[X?,^/jl]^Qi^Pl)6iX?-l) (80) 

for L = 1. We note that the case of 0^ = 1 is a special case where no disorder is 
present and Tq^i = . . . = Tl^i^l since Q{ipc) = ^{i^c ~ ii^c)) for all k. Thus for 
(pd = 1, > 1 is sufficient to show that IRSB solution is unstable at x. To 

evaluate (A^''(x)) for general values of (pa < 1, solving Eq. fl65|) gives more reliable 
results than population dynamics. This is because for large L and nonvanishing values 
of 0s, such as those close to 0^""", the presence of the factor H^i ^A;,a in Xf^{x) requires 
an extremely large population for a finite fraction of nonzero Xf^{x) in the pool of 
population dynamics. 

We show {X^l\x)) as a function of L in Fig. fTBT a) for 0^ = 0.9 from X7 = to 
xj = 4.5, and define x^^^ by 

lim ^Mp-H = 1. (81) 

From Fig. [T6]^a), the IRSB solution for 0^ = 0.9 is stable against type I instability when 
x'y < x'-^^'y f» 3.07. 



6.4-2. Type II instability: Fragmentation of states In the IRSB phase, vertex k is 
characterized by a cavity energy function in each single state. In the 2RSB 
phase, some IRSB single states split into different states in which some vertices are 
characterized by more than one cavity energy functions. In the example illustrated in 
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Fig. [13 a state with = is fragmented to states with E^ = E^ and E^ = E^ 
respectively on entering the 2RSB phase. 

As the cavity energy functions of some vertices are modified during fragmentation, 
we examine the probabihty of changes in E^ of a node i due to changes in Ej among 
its descendent nodes j, the so-called proliferation of bugs [36]. In the singlet regime, we 
denote tt^I^^ as the joint probability that vertex — )■ j is in the S state in the absence of 
bugs, and in the C state in the presence of a small number of bugs. From the recursion 
Eq. ( 173|) . we note that contributions to vr^I^* come from the case that all descendents of 
i are in the S state , and one of them changes to C state in the presence of bugs. Hence 

k&Cj\{i} leCj\{i,k} 

Similarly, contributions to irl^l come from the case that only one descendent of j is in 
the C state which changes to S state in the presence of bugs. Hence 

e 

keCj\{i} ieCj\{i,k} 
We define the matrix Vk-^j{x) to be 



-xu^/2 



4) 



f Yl ^^^^e"^ ^ 

leCj\{i,k} 

n 

V iec,\{i,k} ) 

The instability of the IRSB solution against fragmentation of states can be thus 
examined by considering the maximum eigenvalue of the products of matrices V . We 
measure the quantity 

\f\x) = {K- lfX{VL-l^LVL-2^L^l ■ ■ ■ (85) 

averaged over quenched disorders, where X(. . .) is defined as the maximum eigenvalue of 
the matrix. The IRSB solution is stable against the type II instability if {\^J^\x)) — 
as L — )■ oo. {\^J^\x)) can be solved by population dynamics of Eq. ( I73|) or by solving 
for the distribution on Vl[^^l\'^c] analogous to Eq. ( |79l) . 

We show {\^^\x)) as a function of L in Fig. flGT b) for 0^ = 0.9 from X7 = to 
X7 = 4.5, and define x^^^'^ similarly as Eq. ( IHTj) . From Fig. HGT b). the IRSB solution for 
0rf = 0.9 is stable against type II instability when X7 < x'-^^^7 ~ 2.93. 

x^^'^ and x*^^^-* from 0^ = 0.78 to 1, together with Xs and Xd obtained from the 
complexity curve, are plotted in Fig. flGT c) as a function of (f)d- Reliable results for 
(j)d ~ (p^^ = 0.75 are difficult to obtain as the continuous component in Q{ipc) becomes 
extremely small (see for instance. Fig. [T2|) and extremely high precision is required. 
For (pd lower than ^ 0.88, the IRSB ansatz is stable in the range of x^^^^ < x < x^^\ 
corresponding to the shaded region in Fig. [T6t^c). Both Xg and Xd are found below 
the shaded region, indicating the instability of the physical segment of S. The stable 
range of x lies on the negative segment of S, which implies that states with vanishing 
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Figure 16. (a) {X^p (x)) and (b) (A^^^(a;)) as a function of L for systems with K = 3 
and (j)d — 0.9, from ^7 = to 0:7 = 4.5. (c) a;^^^7, a;(^^'7, Xgj and Xd'y as a function of 
(pd for K — 3. The shaded region corresponds to the stable range of the IRSB ansatz. 

probability are stable in the IRSB ansatz. For (pd higher than ^ 0.88, the shaded region 
disappears and all complexity curves are unstable. We thus conclude that the IRSB 
quantities evaluated at Xs and Xd are unstable. Nevertheless, agreement with simulation 
results show that the IRSB ansatz is a good approximation of the physical picture of 
the system. 

7. Conclusion 

In this paper, we have studied the source location problem on transportation networks. 
As the formulation involves continuous variables, the cavity fields are represented 
by the cavity energy functions which satisfy the piecewise quadratic ansatz. The 
ansatz decomposes the cavity energy functions into composite functions, and effectively 
parametrizes them by the energy minimum of each composite function. 

This enables us to obtain a small closed set of cavity energy functions in the singlet 
regime, which greatly simplifies the functional RS recursions to simple recursions of 
probabilities. Physical results such as the average energy and the fraction of soft nodes 
are obtained and have excellent agreement with simulations when the fraction of deficient 
nodes is small. We examined the stability of the RS solution and derive the AT-line 
for the transition to the RSB phase. In the high connectivity limit, such results are 
consistent with the RS instability obtained in the vertex cover problem for Poissonian 
graphs. 

Though the solution of the IRSB ansatz corresponds to a stable distribution of 
functionals, which is in general infeasible to solve, the closed sets of cavity energy 
functions greatly simplify the IRSB recursions which make the IRSB solution feasible. 
We remark that the analysis is applicable to regimes other than the singlet regime of 
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the system, such as the commensurate points of the doublet regime being considered in 



Appendix A In other cases, closed sets with a large number of functions are found and 
the IRSB solution may once again become computational infeasible. We expect that 
the present techniques are applicable to other problems where closed sets of cavity fields 
exist. 

Appendix A. The doublet Regime 

Appendix A.l. The Closed Set of Cavity Energy Functions and the Simplified RS 
Recursion 

Apart from the closed set of cavity energy functions E^ in the singlet regime, a closed 
set with countable elements of E^ can also be found in other regimes. In the doublet 
regime, we find that a small number of E^ spans a closed set at the commensurate 
points 

7 = mn, (A.l) 

with the rational number m > 2. These commensurate points correspond to the values 
of 



(A.2) 

(m — 1)K + m + 1 

at which abrupt jumps in the fraction of source nodes are found in the range a/21/25 < 
< \f^ij2 in Fig. [2] for K = ?>. Each value of rational m corresponds to a switch of 

energetic stability from one configuration of consumer nodes to another. 

To find the closed set of E^ in the doublet regime, we only have to consider the 

composite functions /o, f\ and /2. E\{yi^ is thus denoted as 

EX[y^) = {dld\,d\). (A.3) 

All other composite functions fny, with > 3 have > (i^. and their corresponding 
configurations are not stable in the doublet regime. We thus consider only the recursion 
A^(0, . . . , 0) = 1 and A^(l, 0, . . . , 0) = 2. As an illustration, we consider the case m = 2 
and 0rf = 1 where closed set of E^ is shown in Fig. lAll The cavity energy change AE^ 
from Eq. fl30|) . the constant terms (ig, d\ and d\ from Eqs. (I22l) and fl29l) are simplified 
to 

k&Cj\{i} 

di = 0, (A.4) 
d{ =-2k- ^ mm.{Q, d\,dl), 

k&Cj\{i} 



di =~K+ min d'l - ^ min(0, dj', ^2), 



k&Cj\{i} 
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Figure Al. A closed set of cavity energy functions Ek{y) &i K = i, = 1 in the 
doublet regime. The forms of corresponds to (a) the a state, (b-c) the c state, (d) 
the d state and (e-h) the s state. 



where the zero-point of E is made at the minimum values of /o for convenience of 
analysis. The energy change of AEj of adding vertex j is obtained by comparing the 
energies of the S states of the vertex and its descendents, and can be shown to produce 
the same physical results as the full energetic comparison. From the recursion relations, 
the full closed set of is found to be 

^<r(?/) = (0,-2«:,(-l + r)/t) 
Er(y) = (0,-/s:,(-l + r)/^) 
Ej(y) = (0,r/t,-l/s:) 

E'^{y) = {Q,qK,rK) (A.5) 

with integers g, r > 0. E"^ corresponds to cavity states with a strong preference to be 
singly consuming, i?^, E^^ and E^ correspond to cavity states which are respectively 
singly consuming, doubly consuming and resource providing. They are denoted as the 
a, c, d and s states. Note that the integer r > 1 in Eq. ( ]A.5I) may correspond to the 
form of E^ where /i and /2 are not relevant (see for instance Fig. lAll (f) and (g)). For 
= 3, there are two forms of E^ with r = and r > 1 and four forms of E^ with (g, r) 
= (> 1, 1), (> 1,0), (0, > 1) and (0,0). Hence the closed set of E^ has eight forms of 
E^ as shown in Fig. lAll 

Next, we consider the disordered case 0^ < 1. The recursion relations can be 
simplified in terms of ipa-, '^c, '4'd and ips, corresponding to the probabilities for a vertex 
to be in the a, c, d and s states, as given by 

keCj\{i} 
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Figure A2. (a) The K-(f)d phase diagram and the AT line (dashed line with circles) 
in the doublet regime with li-i y^KiK ~ 1)/{K + 3). The dotted line is the AT line 
in the singlet regime, (b) The complexity function S(e) from the IRSB solution for 
K = 3 with u~-^ = 1. e' = e/(w^/2). Inset: The negative segment of the complexity 
function shown in the limit x — > (X). Each symbol on the line indicates an increase of 
1.2 in XK. 



k€Cj\{i} l&Cj\{i,k} 
k€Cj\{i} l<^Cj\{i,k} 

and V^r' = 1 -^r' - V^r' - V^d^' In the RS phase, we set ipi-^', -^y', ipl-^' = 0, 1. 
The average optimized energy in the RS ansatz is obtained by evaluating the full energy 
change A^node and Af^unk resulting from the addition of new nodes and new links. The 
expression for Af^nodo and Af^unk from Eq. (jH]) can be simplified as 

2 

^^node = y + min [ E t^H^, di, di), -2k, -k + min d{ | 

+ 8k,,A E niin(0, d{, di). 

ASiin^ = min [o, 4, df, _ + 4 + <, 4, df^ . (A.7) 
Appendix A. 2. The Instability of the RS Ansatz 

To obtain the AT line between the RS and the RSB phases, we consider the variations 
6ip'^^^, 6ip'^^^, 5%l)l^^ and 5^1)^^^^ in the recursions of probabilities in Eq. (lA.6p . with 



Sip^~^^ _l_ Sip'^^^ + Sip^~^^ + Sipg^^ = 0. Though the c and d states lead to different 
degeneracies of the full states, the two states play the same role in the recursion 
relations and lead to the same AE, A£^node and Af^unk- Combining ipc^'' and V^^^"' 
in the recursions, the variations 6ip^~^^, 5{ijj^~^^ + V'd^"') and 5ijj^~^^ depends only on 
^ip'^-^j of the descendents. We thus write down the AT condition for the disordered case 
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of the doublet regime, 

= (K-i){K- 2)0,(1 - (a.s) 

where we have apphed the relation of {{'ip^'^-'Y) = (^^^"') derivation. Simple 

algebra leads to the following form of the AT line 

^AT ^ ]_ 

{K - 1){K - 2)(1 - (^.)AT)((V;,)AT)i^-3' 



2K{K-2) ' ^ ' 

where (ips)^"^ is the average value of ips on the AT line. The K-cf)^ phase diagram is 
shown in Fig. IA2( a). which is compared to the phase diagram of the singlet regime in 
Fig. [H] inset. The comparison suggests that the RSB phase shrinks when decreases 
from the singlet regime to the doublet regime. 



Appendix A. 3. The IRSB solution 

To obtain the IRSB solution, we follow the approach adopted in the singlet regime and 
evaluate the partition functions Snode and Sunk using the full energy change £^node and 
^^link- We solve the IRSB solution for the case without disorder (i.e. (pd = 1) with the 
IRSB restriction on delta functional form of Q{P[£'^]}, i.e. Q{4'a,i'c,4'd) = ^("^a — 
{4'a))S{i^c — {i^c))S{i^d — ii^d))- The complexity function S is obtained from the IRSB 
solution for K = 3, and is shown in in Fig. IA2( b). Compared with the singlet regime, S 
in both cases have similar form and similar maximum values. The physical segment of 
the curve is shown by the solid segment. Note that e= (1 — /c — /d)/2 + /c/6 + /x)/3 and 
the fraction fc and ft> of singly and doubly consuming nodes are not uniquely determined 
by the energy e. The negative segment of the complexity function approaches the limit 
of e = 1/3 as X — )■ oo, as shown in the inset of Fig. IA2( b). It corresponds to the 
lowest possible energy on graphs with special structures, which occur with vanishing 
probability in the thermodynamic limit, e = 1/3 implies 4/c + 3/© = 2, suggesting 
the modulation limit of /c = 1/2 when /x) = as in the singlet regime, or the limit of 
f-D = 2/3 when /c = where two-third of nodes are doubly consuming. 
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